Classification of Variants across paired longitudinal samples in the PBTA Cohort

Set up

suppressPackageStartupMessages({
  library(tidyverse)
})

Directories and File Inputs/Outputs

# Detect the ".git" folder -- this will be in the project root directory
# Use this as the root directory to ensure proper sourcing of functions 
# no matter where this is called from
root_dir <- rprojroot::find_root(rprojroot::has_dir(".git"))
analysis_dir <- file.path(root_dir, "analyses", "tmb-vaf-longitudinal")
results_dir <- file.path(analysis_dir, "results")
input_dir <- file.path(analysis_dir, "input")
files_dir <- file.path(root_dir, "analyses", "sample-distribution-analysis", "results")

# Input files
genomic_paired_file <- file.path(files_dir, "list-genomic-paired.tsv") 
tmb_vaf_file <- file.path(results_dir, "tmb_vaf_genomic.tsv")
palette_file <- file.path(root_dir, "figures", "palettes", "oncoprint_color_palette.tsv")

# File path to plot directory
plots_dir <-
  file.path(analysis_dir, "plots")
if (!dir.exists(plots_dir)) {
  dir.create(plots_dir)
}

source(paste0(root_dir, "/figures/scripts/theme.R"))

Read in data and process

tmb_vaf_df <- readr::read_tsv(tmb_vaf_file, guess_max = 100000, show_col_types = FALSE) %>% 
  filter(!tmb >= 10) %>% 
  select(Kids_First_Biospecimen_ID, Variant_Classification, gene_protein, mutation_count,   region_size, tmb, VAF)

genomic_paired_df <- readr::read_tsv(genomic_paired_file, guess_max = 100000, show_col_types = FALSE) %>% 
  left_join(tmb_vaf_df, by = c("Kids_First_Biospecimen_ID"))

# Which patient samples don't have TMB?
#genomic_paired_df %>% 
#  filter(is.na(tmb)) %>% 
#  unique() %>% 
#  regulartable() %>%
#  fontsize(size = 12, part = "all")

# Vector to order timepoints
td_order <- c("Diagnosis", "Progressive", "Recurrence", "Deceased", "Second Malignancy", "Unavailable")

df <- genomic_paired_df %>% 
  filter(!is.na(tmb)) %>% 
  mutate(tumor_descriptor = factor(tumor_descriptor),
         tumor_descriptor = fct_relevel(tumor_descriptor, td_order))


# Let's count #samples per cancer groups and timepoints 
timepoint_n_df <- df %>% 
  count(cancer_group, tumor_descriptor) %>% 
  dplyr::mutate(tumor_descriptor_cancer_n = glue::glue("{cancer_group}_{tumor_descriptor}  (N={n})")) %>% 
  dplyr::rename(timepoint_n = n) 


# Create df to use for plots
df_plot <- df %>% 
  left_join(timepoint_n_df, by = c("tumor_descriptor", "cancer_group")) #%>%
  #filter(!timepoint_n <= 2) 

Type of non synonymous mutations per cancer group

# Read color palette
palette_df <- readr::read_tsv(palette_file, guess_max = 100000, show_col_types = FALSE) 

# Define and order palette
palette <- palette_df$hex_codes
names(palette) <- palette_df$Variant_Classification

Warning: Unknown or uninitialised column: `Variant_Classification`.

# Define label for plots
Alteration_type <- df_plot$Variant_Classification

# Define ylim
ylim <- max(df_plot$tmb)

What type of alterations we observe per tumor descriptor?

# Create bxp
print(ggpubr::ggboxplot(df_plot, 
                        x = "tumor_descriptor", 
                        y = "tmb", 
                        color = "Variant_Classification",
                        #add = "jitter",
                        palette = palette) +
        theme_Publication() + 
        scale_y_continuous(limits = c(0, ylim)) +
        ylab("TMB") +
        xlab("Timepoint") +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = "Alteration_type_timepoints.pdf", 
       path = plots_dir, 
       width = 15, 
       height = 8, 
       device = "pdf", 
       useDingbats = FALSE)

What type of alterations we observe per tumor descriptor in each cancer group?

# Create bxp
print(ggpubr::ggboxplot(df_plot, 
                        x = "tumor_descriptor", 
                        y = "tmb", 
                        color = "Variant_Classification",
                        #add = "jitter",
                        palette = palette) +
        facet_wrap(~cancer_group) +
        theme_Publication() + 
        ylab("TMB") +
        xlab("Timepoint") +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = "Alteration_type_cancer_group.pdf", 
       path = plots_dir, 
       width = 25, 
       height = 18, 
       device = "pdf", 
       useDingbats = FALSE)

What type of alterations we observe per tumor descriptor in each cancer group (add _n))?

cg <- as.character(unique(df_plot$cancer_group))
cg

 [1] "CNS Embryonal tumor"                       "Embryonal tumor with multilayer rosettes"  "Low-grade glioma"                         
 [4] "Ependymoma"                                "Medulloblastoma"                           "Atypical Teratoid Rhabdoid Tumor"         
 [7] "Diffuse midline glioma"                    "High-grade glioma"                         "Ganglioglioma"                            
[10] NA                                          "Meningioma"                                "Glial-neuronal tumor"                     
[13] "Pilocytic astrocytoma"                     "Neuroblastoma"                             "Ganglioneuroblastoma"                     
[16] "Schwannoma"                                "Chordoma"                                  "Fibromyxoid lesion"                       
[19] "Neurofibroma/Plexiform"                    "Malignant peripheral nerve sheath tumor"   "Choroid plexus carcinoma"                 
[22] "Adamantinomatous Craniopharyngioma"        "Dysembryoplastic neuroepithelial tumor"    "Ewing sarcoma"                            
[25] "Rosai-Dorfman disease"                     "Hemangioblastoma"                          "Craniopharyngioma"                        
[28] "Diffuse leptomeningeal glioneuronal tumor"

# Loop through variable
for (i in seq_along(cg)){
  print(i)
  df_sub <- df_plot %>%
      filter(cancer_group == cg[i])
  
  # Create bxp
  print(ggpubr::ggboxplot(df_sub, 
                        x = "tumor_descriptor_cancer_n", 
                        y = "tmb", 
                        color = "Variant_Classification",
                        palette = palette) +
        theme_Publication() + 
        ylab("TMB") +
        xlab("Timepoint") +
        labs(title = paste(cg)) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)))

  # Save the plot
  #name <- paste0(cg[i], "-Alteration_type_cancer_group.pdf")
  #ggsave(filename = name, 
  #       path = plots_dir, 
  #       width = 12, 
  #       height = 8, 
  #       device = "pdf", 
  #       useDingbats = FALSE)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28

What type of alterations we observe per tumor descriptor in each cancer group and timepoint model?

tm <- as.character(unique(df_plot$timepoints_models))
tm

 [1] "Dx-Rec"         "Dx-Pro"         "Dx-Dec"         "Pro-Dec"        "Rec-SM"         "Pro-Rec"        "Dx-SM"          "Pro-Rec-Dec"    "Dx-Pro-Rec"    
[10] "Rec-Dec"        "Dx-Pro-Rec-Dec" "Dx-Rec-Dec"     "Dx-Pro-Dec"     "Dx-Unv"        

# Loop through variable
for (i in seq_along(tm)){
  print(i)
  df_sub <- df_plot %>%
      filter(timepoints_models == tm[i])
  
  # Create bxp
  print(ggpubr::ggboxplot(df_sub, 
                        x = "tumor_descriptor", 
                        y = "tmb", 
                        color = "Variant_Classification",
                        palette = palette) +
        facet_wrap(~cancer_group) +
        theme_Publication() +
        ylab("TMB") +
        xlab("Timepoint") +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))

   # Save the plot
   #name <- paste0(timepoints_models[i], "-Alteration_type_cancer_group.pdf")
   #ggsave(filename = name, 
   #    path = plots_dir, 
   #    width = 12, 
   #    height = 8, 
   #   device = "pdf", 
   #   useDingbats = FALSE)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14

sessionInfo()

R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.5.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggthemes_4.2.4  flextable_0.9.2 lubridate_1.9.2 forcats_1.0.0   stringr_1.5.0   dplyr_1.1.2     purrr_1.0.1     readr_2.1.4     tidyr_1.3.0    
[10] tibble_3.2.1    ggplot2_3.4.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] fontquiver_0.2.1        bit64_4.0.5             rprojroot_2.0.3         tools_4.2.3             backports_1.4.1         bslib_0.5.0            
 [7] utf8_1.2.3              R6_2.5.1                colorspace_2.1-0        withr_2.5.0             tidyselect_1.2.0        bit_4.0.5              
[13] curl_5.0.2              compiler_4.2.3          textshaping_0.3.6       cli_3.6.1               xml2_1.3.5              officer_0.6.2          
[19] fontBitstreamVera_0.1.1 labeling_0.4.2          sass_0.4.7              scales_1.2.1            askpass_1.2.0           systemfonts_1.0.4      
[25] digest_0.6.33           rmarkdown_2.23          gfonts_0.2.0            pkgconfig_2.0.3         htmltools_0.5.5         fastmap_1.1.1          
[31] rlang_1.1.1             rstudioapi_0.15.0       httpcode_0.3.0          shiny_1.7.4.1           jquerylib_0.1.4         farver_2.1.1           
[37] generics_0.1.3          jsonlite_1.8.7          vroom_1.6.3             zip_2.3.0               car_3.1-2               magrittr_2.0.3         
[43] Rcpp_1.0.11             munsell_0.5.0           fansi_1.0.4             abind_1.4-5             gdtools_0.3.3           lifecycle_1.0.3        
[49] stringi_1.7.12          yaml_2.3.7              carData_3.0-5           plyr_1.8.8              parallel_4.2.3          promises_1.2.0.1       
[55] crayon_1.5.2            hms_1.1.3               knitr_1.43              pillar_1.9.0            ggpubr_0.6.0            uuid_1.1-0             
[61] ggsignif_0.6.4          reshape2_1.4.4          pkgload_1.3.2.1         crul_1.4.0              glue_1.6.2              evaluate_0.21          
[67] fontLiberation_0.1.0    data.table_1.14.8       vctrs_0.6.3             tzdb_0.4.0              httpuv_1.6.11           gtable_0.3.3           
[73] openssl_2.1.0           cachem_1.0.8            xfun_0.39               mime_0.12               xtable_1.8-4            broom_1.0.5            
[79] rstatix_0.7.2           later_1.3.1             ragg_1.2.5              timechange_0.2.0        ellipsis_0.3.2